R Notebook

This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.

Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.

library(qiime2R)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(tidyr)
library(microbiome)

Loading required package: phyloseq
Loading required package: ggplot2

microbiome R package (microbiome.github.com)
    


 Copyright (C) 2011-2019 Leo Lahti, 
    Sudarshan Shetty et al. <microbiome.github.io>


Attaching package: ‘microbiome’

The following object is masked from ‘package:ggplot2’:

    alpha

The following object is masked from ‘package:base’:

    transform

library(compositions)

Welcome to compositions, a package for compositional data analysis.
Find an intro with "? compositions"


Attaching package: ‘compositions’

The following objects are masked from ‘package:stats’:

    anova, cor, cov, dist, var

The following objects are masked from ‘package:base’:

    %*%, norm, scale, scale.default

library(Dict)
library(reshape2)

Attaching package: ‘reshape2’

The following object is masked from ‘package:tidyr’:

    smiths

library(heatmaply)

Loading required package: plotly

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout

Loading required package: viridis
Loading required package: viridisLite
Registered S3 method overwritten by 'dendextend':
  method     from 
  rev.hclust vegan
Registered S3 method overwritten by 'seriation':
  method         from 
  reorder.hclust vegan

======================
Welcome to heatmaply version 1.3.0

Type citation('heatmaply') for how to cite the package.
Type ?heatmaply for the main documentation.

The github page is: https://github.com/talgalili/heatmaply/
Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
You may ask questions at stackoverflow, use the r and heatmaply tags: 
     https://stackoverflow.com/questions/tagged/heatmaply
======================


Attaching package: ‘heatmaply’

The following object is masked from ‘package:compositions’:

    normalize

library(KEGGREST)

#This code block is to stage all the files needed for visualization


#Stage taxonomy file using Qiime2R
taxa_vsearch <- parse_taxonomy(taxonomy = read_qza('../taxonomy_vsearch/classification.qza')$data)
taxa_blast <- parse_taxonomy(taxonomy = read_qza('../taxonomy_blast/classification-blast.qza')$data)

#Stage metadata file
sample_metadata <- read.csv('../read-files/metadata.tsv', sep='\t')

#Stage unstratified/strat kegg pathway files
unstrat_kegg <- read.csv('../picrust2-strat/out/KEGG_pathways_out/path_abun_unstrat.tsv', sep = '\t')
strat_kegg <- read.csv('../picrust2-strat/out/KEGG_pathways_out/path_abun_strat.tsv', sep = '\t')

#Stage unstrat/strat KEGG ortholog files 
ko_unstrat <- read.csv('../picrust2-strat/out/KO_metagenome_out/unstrat_ko_metagenome.tsv', sep = '\t')
ko_strat <- read.csv('../picrust2-strat/out/KO_metagenome_out/strat_ko_metagenome.tsv', sep = '\t')

#Open selected KOs (for plotting relevant Kegg pathways)
deg_kos <- read.csv('./select_kos.txt', sep = '\t', header = TRUE)
rownames(deg_kos) <- deg_kos[,1]
deg_kos[,1] <- NULL

#open the pathway to ko file
pth_2_ko <- read.csv('~/bioinfo_pipelines/picrust2-2.4.2/picrust2/default_files/pathway_mapfiles/KEGG_pathways_to_KO.tsv', header =FALSE, sep='\t')
#Filter to only Select KOs

pth_2_ko <- pth_2_ko[pth_2_ko$V1 %in% rownames(deg_kos),]

#Retrieve Sample ids so I can call it as a variable 
sample_ids <- c(as.character(sample_metadata$Sample.id))

#retrieve BRITE hierarchy from KEGG using KEGGREST

keggGet(dbentries = 'br08901')

Error in .getUrl(url, .flatFileParser) : Not Found (HTTP 404).

#This code block is to load the required files to plot NSTI vs relative abundance
abs_count_otu <-  read_qza(file = '../dada2/table.qza')$data
rel_count_otu <- as.data.frame(apply(abs_count_otu,2,function(x){x/sum(x)}))
nsti_asvs <- as.data.frame(read.csv('../picrust2-strat/out/marker_predicted_and_nsti.tsv.gz', sep='\t', row.names = 1))


rel_nsti <- merge(nsti_asvs,rel_count_otu,by="row.names",all.x=TRUE)


rel_nsti <- melt(rel_nsti, measure.vars = sample_ids, variable.name = 'sample.id', value.name = 'rel.abundance')


#Group accding to consortia



ggplot_rel_plot <- ggplot(data=rel_nsti,mapping= aes(x=metadata_NSTI, y=rel.abundance, color=sample.id)) + geom_point()

plot(ggplot_rel_plot)